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Abstract. We study a two-dimensional two-component Fermi gas with 
attractive or repulsive short-range interactions at zero temperature. We 
use Diffusion Monte Carlo with Fixed Node approximation in order to 
calculate the energy per particle and the opposite-spin pair distribution 
functions. We show the relevance of beyond mean field effects and verify 
the consistency of our approach by using Tan's Contact relations. 



1 Introduction 

Low dimensional configurations of degenerate Fermi and Bose gases are being object 
of many experimental and theoretical studies [T] as a means of simulating strongly 
correlated systems and as interesting systems per se, from a fundamental point of 
view. 

Many important areas of investigation on Fermi gases are being extended to 
two-dimensional (2D) configurations, which have already been pursued in the three- 
dimensional (3D) case, such as the BCS-BEC crossover in a superfluid gas with reso- 
nantly enhanced interactions and the possible onset of itinerant ferromagnetism in a 
gas with repulsive interactions. Other interesting phenomena that should be promi- 
nent in the 2D case are the Fulde-Ferrell-Larkin-Ovchinnikov superfluid state, the 
ultracold gases analogue of the quantum Hall effect in presence of non-abelian gauge 
fields and the long range correlations induced by dipolar interactions. 

In particular, 2D ultracold Fermi gases have received a lot of attention in recent 
years, having been realized using highly anisotropic pancake-shaped potentials. Den- 
sity profiles of the clouds have been measured using in situ imaging [213] : the single 
particle spectral function has been measured by means of rf spectroscopy [I], the 
presence of many-body polaronic states and their transition to molecular states have 
been characterized [5]. 

On the theoretical side, the solution of the BCS equations for the 2D attractive 
Fermi gas has been first investigated in [6] and later in [7]. The perturbative anal- 
ysis of the 2D repulsive Fermi gas has been performed in [17]. Recent experiments 
are in a regime where beyond mean-field contributions become relevant. In such a 
respect, many results are now available especially for the highly imbalanced case [5], 
both for the ground-state and for the so-called Upper Branch (UB). The 2D case is 
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particularly challenging from the point of view of theory, being a marginal case in 
field theories, when the leading dependence of the observables on the coupling is non 
algebraic. Recently we have obtained the first determination using Quantum Monte 
Carlo methods of the equation of state at T = of a balanced homogeneous 2D Fermi 
gas in the BCS-BEC crossover [5]. Such an equation of state has been positively com- 
pared to experiments using the local density approximation [3J, it has been used to 
discuss recent results on collective modes in a 2D Fermi gas pPQJ and to judge the 
temperature dependence of the Contact in 2D trapped gases [11) . 

In this article we report new results concerning the equation of state of the repul- 
sive gas and its density-density correlation function. Moreover we emphasize the need 
of a consistent choice for the coupling constants and the coefficients of the beyond 
mean-field contributions to energy, both in the strongly interacting molecular regime 
and in the weakly interacting Fermi liquid regime. In Section [2] we introduce the 
model potentials and we discuss the trial nodal surfaces used in the Quantum Monte 
Carlo method; in Section [3J we present the equation of state of the weakly attractive 
or repulsive gases, we show results for the Upper Branch of the attractive gas and 
we discuss the equation of state of the composite bosons in the molecular regime; 
finally in Section 2] we discuss the extraction of the Contact from the density-density 
correlation function and we show results for the Contact of the repulsive gas. 



We consider a homogeneous two-component Fermi gas in 2D described by the Hamil- 
tonian 



where m denotes the mass of the particles, ... and label, respectively, 

spin- up and spin-down particles and Nf = iVj. = N/2, N being the total number of 
atoms. We model the interspecies interatomic interactions using three different types 
of model potentials: an attractive square- well (SW) potential V(r) = — Vq for r < Rq 
(Vq > 0), and V(r) — otherwise; a repulsive soft-disk (SD) potential V(r) = V\ for 
r < R\ (Vi > 0), and V(r) = otherwise; and a hard-disk (HD) potential V(r) = oo 
for r < i?2 and V(r) = otherwise. Due to the logarithmic dependence on energy of 
the phase shifts in 2D, different definitions of the scattering length have been used [12] : 
we define the scattering length a 2 o — R2 in the case of the HD potential, so that 
for the SD potential one gets a 2 D — Ri e~ Io ^ Kl ^ Kl11 ^ and for the SW potential 
aw = Ro e J °( K °>i K ° Jl ( Ko \ where J n) and I m are the Bessel and modified Bessel 
functions of first kind and ko = \J V^mR 2 / h 2 , K\ = \JVimR\/h 2 . In order to ensure 
the diluteness of the attractive gas we use uRq = 10~ 6 , where n is the gas number 
density. The Fermi wave vector is defined as kp — s/2~7vn, and provides the energy 
scale Ef — h 2 kp/2m. In order to assure universality in the repulsive case we check 
that the results for the HD potential and for the SD potential with different values 
of i?i are compatible. 

For the attractive SW potential in 2D the scattering length is always non negative 
and diverges at null depth and at the zeros of J\ , corresponding to the appearance of 
new two-body bound states in the well. For a SW potential therefore a bound state 
is always present, no matter how small the attraction is. The shallow dimers have 
size of order a 2 D and their binding energy is given by Eh — — Ah 2 / {ma\ D e 21 ), where 
7 ~ 0.577 is Euler-Mascheroni's constant. A different and widely used definition of 
the 2D scattering length is b = a2De 7 /2, such that £& = —H 2 /mb 2 , analogously to 
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Fig. 1. (a) Scattering length of the SW potential, as a function of the depth Vo/ef, in the 
regime used in the simulations. The dashed line corresponds to \jkFaio = 1. (b) Scattering 
length of the SD potential, as a function of the barrier V\. The dotted line corresponds to 
one of the values of R\ used in the simulations. 



the 3D case. The dependence of a,2D on the depth Vp in the region where the well 



supports only one bound state is shown in Fig. 1(a) the dashed line indicates the 
value of scattering length at which the bare molecules have a size comparable to the 
mean interparticle distance 1/kp. The region /cf02D 2> 1 corresponds to the BCS 
regime where interactions are weak and dimers are large and weakly bound, while 
kFa>2D <C 1 corresponds to the BEC regime of tightly bound composite bosons. The 
regime in which the scattering length diverges and a bound state appears (unitarity 
limit) is trivial in 2D because it corresponds to the non-interacting case; instead the 
resonant regime corresponds to the region kpa2D ~ 1; this can be seen at the two-body 
level by considering the low-energy scattering amplitude f(k) = 27r/[log(2/fca2De 7 ) + 
in/ 2] [T2], which is enhanced for k ~ I/oqd (with logarithmic accuracy). 

For the repulsive SD potential the scattering length is always positive and smaller 
than Ri, like in 3D (see Fig. 1(b) ). Although the real physical potentials always have 
an attractive part, the purely repulsive models that we use can be useful in describing 
the cases in which the scattering length is small and positive compared to the mean 
interparticle distance, when it is possible to prepare metastable repulsive gas-like 
states without a significant production of molecules [5] . 



2.1 Monte Carlo simulations 

We use the fixed-node diffusion Monte Carlo (FN-DMC) method. This numerical 
technique solves the many-body Schroedinger equation by an imaginary time pro- 
jection of an initial guess of the wavefunction. Provided that the initial guess has a 
finite overlap with the true ground-state, this method provides the exact energy of a 
systems of bosons, with a well controllable statistical error. For fermions, FN-DMC 
yields an upper bound for the ground-state energy of the gas, resulting from an ansatz 
for the nodal surface of the many-body wave function that is kept fixed during the 
calculation (see Refs. [14]). 

The fixed-node condition is enforced using an initial and guiding trial function 
that we choose of the standard form V't(R) = ^s(R)^(R), namely the product 
of a purely symmetric and a purely antisymmetric term. $>a satisfies the fermionic 
antisymmetry condition and determines the nodal surface oftpx, while <Ps is a positive 
function of the particle coordinates and is symmetric in the exchange of particles with 
equal spin (Jastrow function). 

Two opposite regimes are described by the <£>a component. The deep attractive 
BEC regime, where the opposite-spin fermions are expected to pair into a condensate 
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of dimers, can be described by an antisymmetrized product #^(R) = A (<t){rii')4>(r22>)---<t l (' r N r N i )) 
of pairwise orbitals <j) corresponding to the two-body bound state of the potential 
V(r). This wavefunction has been proposed by Leggett as a projection of the grand- 
canonical BCS function to a state with a finite number of particles, later extended to 
the polarized case [13] and extensively used in 3D Quantum Monte Carlo simulations 
[14) . The weakly interacting regime, where a Fermi liquid description is expected to 
be valid, can be instead described by a typical Jastrow-Slater (JS) function with 
<Pa(R-) = Df(Nf)Di(Ni), namely the product of the plane-wave Slater determinants 
for spin-up and spin-down particles. This description is expected to hold both in the 
weakly repulsive branch and in the attractive BCS regime of a weakly interacting gas 
where the effect of pairing on the ground-state energy is negligible. 

The symmetric part is chosen of the Jastrow form ^s(R) = /^(R) J-j-)-(R) J^(R). 
The diluteness of the gas allows us to consider just two-body Jastrow functions, with 

■%(R) = n M '/n(»v), ^tt(R) = Ui^Mnj), J u pQ = Il v<j > fufoj'), where 

two-body correlation functions of the interparticle distance have been introduced, 
which aim at reducing the statistical noise by fulfilling the cusp conditions, namely the 
exact behavior that is expected when two particles come close together. In particular 
in the weakly interacting attractive or repulsive regimes we set fn = = 1, while 
fol is set equal to the analytical ground-state solution of the two-body Schroedinger 
equation in the center-of-mass frame, with the bare two-body potentials. In the 
strongly attractive regime instead ffi = 1, since the opposite-spin short-range corre- 
lation is already accounted for in the BCS orbitals; for the parallel spin correlations 
we use the ground-state solution of an effective two-body problem, consisting of a SD 
interaction with scattering length a e // = 0.6a2_o, motivated by the expected interac- 
tion between molecules [12]. This choice greatly reduces the variance of the sampled 
energy. 

In this paper we also study the Upper Branch in the strongly attractive regime, 
which corresponds to a (metastable) gas of repulsive fermions; being a many-body 
excited state, we only perform a Variational Monte Carlo (VMC) simulation, without 
imaginary time projection. In this case the wavefunction is taken of the JS type, 
with the ff± functions equal to the first excited state of the two-body SW problem, 
aiming at enforcing orthogonality to the molecular ground-state. The nodal surface in 
this case arises also from the "Jastrow" function. This unusual non positive definite 
factor has already been used in the context of 3D Fermi gases for the investigation of 
itinerant ferromagnetism [15) . 

Simulations are carried out in a square box of area L 2 = N/n with periodic 
boundary conditions (PBC). In order to fulfill those, all radially symmetric two-body 
functions have zero derivative at r = L/2, both in the Jastrow factors and in the 
molecular orbitals in the BCS wavefunction; this is obtained by smoothly matching 
each of the previously described functions to a sum of exponentials f e (r) = c\ + 
C2[exp(— fi(L — r))+exp(— fir)] at some healing length R, with [i and R parameters to 
be optimized. The PBC also select a specific set of compatible plane- waves to be used 
in the JS wavefunction. It is known that finite-size effects of the JS wavefunction can 
be strongly depressed by considering numbers of fermions which provide a maximally 
symmetric Fermi surface (closed shells): we choose N/2 — 13 and N/2 — 49. Moreover 
we add a correction which can be justified with the theory of Fermi liquids (see [16]). 
namely the energy difference between the finite and infinite system is assumed to be 
the same as for the noninteracting case. Since this cannot assure the elimination of all 
finite-size effects (and in particular we do not consider the role of the effective mass), 
we use this correction also to assess the error-bars, on top of the statistical error. No 
significant finite-size effect is seen when using the BCS trial function. 
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Fig. 2. Energy per particle of a weakly interacting 2D Fermi gas (inset), with the mean- 
field contribution subtracted (main figure). The attractive regime corresponds to g < 0, the 
repulsive regime to g > 0. The results from the simulations with the SW and HD potentials 
are compared to the Fermi liquid theory result @. The dashed line indicates the quadratic 
fit to the SW JS-DMC data. Empty symbols refer to calculations with the SD potential, 
with Ri/a,2D = 2, 4, only shown for two representative values of g. 



3 Energy 

3.1 Weakly interacting regime 

The Fermi Liquid Theory (FL) prediction for the zero temperature equation of state 
of a weakly short-range interacting 2D gas [T7] is the following: 

E FL /N = E FG [l + 2. 9 + (3-41o g 2). 9 2 ] , (2) 

where Efg = h 2 k 2 F /Am = Ef/2 is the energy per particle of the noninteracting 
gas and g is the coupling constant. The peculiar logarithmic dependence of the 
2D scattering amplitude on the available kinetic energy leads to an arbitrary de- 
pendence of the 2D coupling constant g = 1/ \og(EA/Ex) on a reference energy 
Ek, the Ea parameter being relative to the specific choice of the potential. When 
finite-range effects are negligible Ea = 4ft 2 / \ma\ D e 21 ), which is equal to |ej,| in case 
of an attractive potential. Since we consider a weakly interacting Fermi liquid, we 
make the appropriate choice Ek = 2ef. Therefore we use the coupling constant 
g = — 1/ \og(na^ D co) — — l/21og(fc^6) with cq = ire 2 ~ / /2. An important remark is in 
order: while the coefficient of the first order term is independent of Cq, the reported 
value of the coefficient of g 2 depends on the chosen cq, since a modification of cq gives 
rise to higher order contributions to the equation of state. 

In Fig. [5] we show the FN-DMC results with the JS nodal surface. Negative cou- 
plings g < correspond to the ground-state in the weakly attractive regime, while 
positive couplings correspond to either the ground-state of the HD potential or to the 
Upper Branch of the attractive system. In the inset we present the energy, while in the 
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main figure we show the energy with the mean-held result Efg(1 + 2g) subtracted, 
in order to emphasize beyond mean-field contributions. Although the three series of 
data obviously confirm the mean-field prediction, interesting different behaviors are 
found when checking the residual contributions to energy. In the g < regime we fit 
the coefficient A2 of the second order term obtaining A2 = 0.69(4), to be compared 
to 3 - 4 log 2 = 0.227. We find A 2 « 3 - 4 log 2 only if we set c = 27r (as has been 
done in [8]). On the opposite side g > 0, the Upper Branch (UB-VMC) results are 
approximately close to the second order fit done for the g < data, at least for small 
g. The discrepancy is due to the lack of optimization of the UB wavefunction, which is 
problematic since there is not a variational principle for this state. The JS results for 
the SD potential with R\ = 2a 2 D are compatible with those with the HD potential. 
We also perform simulations with shallower SD potentials (R\ — 4o2d), which result 
in departures from universality at g 0.25. Agreement between eq. |2j and the HD 
JS-DMC results is found up to g « 0.1, which corresponds to na 2 ~ 10~ 3 . The be- 
yond mean-field contributions in this paramagnetic phase soon start to decrease. In a 
related system a ferromagnetic transition has been excluded [18] ; further investigation 
concerning the strongly repulsive regime is however beyond the scope of this paper. 
The intriguing difference in the beyond mean-field terms, between the attractive case 
and the purely repulsive case, would also deserve further investigation, hopefully with 
an analytical treatment. 



3.2 BCS-BEC crossover 

We now pass to discuss the main result published in [8] , namely the characterization of 
the composite bosons equation of state in the 2D BCS-BEC crossover. The 2D mean- 
field BCS equations can be analytically solved |6I7] along the BCS-BEC crossover, in 
terms of the interaction coefficient x — |£f,|/2£p. For the BCS order parameter one 
obtains A — 2sf\/x, while for the energy per particle one obtains E/N = Efg0- — 2%) 
and for the chemical potential fi = sp(l — x). For small binding energy one recovers 
the non interacting limit; when x ~ 1, the chemical potential becomes zero and 
then negative, so that the role of the dimers becomes more important; for very strong 
binding the chemical potential of the fermions is equivalent to half the binding energy 
of a molecule, so the system is made of non interacting bosons. Although very useful 
for providing a global self-consistent picture and for setting a stringent variational 
upper bound to the energy per particle, the BCS solution fails in various aspects. In 
the BCS regime it neglects the Hartree-Fock contributions to energy (discussed in the 
previous Section), which are dominant, since the gap is small. In the BEC regime it 
misses the correct interaction energy between the bosons. In general it is not able to 
reproduce the logarithmic dependence of energy on the density, which is typical of 
2D. 

In Fig. [3] we show the FN-DMC results (from [8]) for the equation of state of a 
short-range attractively interacting 2D g function of the interaction parame- 

ter rj = \og(kpa2D) in units of Efg, with Sb/2 subtracted. The BCS wavefunction 
provides a lower energy for values of the interaction parameter rj < 1, while the JS 
function is more favorable for larger values of ij. In the deep BEC regime, besides the 
molecular contribution, the remaining fraction of energy corresponds to the interac- 
tion energy of the bosonic dimers. 

In the BEC regime the FN-DMC results are fitted with the equation of state of 
a 2D gas of composite bosons. The functional form is derived from the analysis of 
Bcane [19] , in the framework of quantum field theory; it also corresponds to previous 
analytical and DMC studies (see Q]5] and references therein). Following Beane we 
introduce the running coupling <?c/(A) in terms of the scattering length of the bosons 
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Fig. 3. Energy per particle of an attractive 2D Fermi gas along the BCS-BEC crossover 
with the bare binding energy of the molecules subtracted. The BEC regime corresponds to 
g — ¥ —oo, the BCS regime corresponds to g — > +oo, while the resonant regime corresponds 
to g ~ 1. The various curves are described in the text. Reprint of the data in [5]. 



aa and the particle density of the bosons nd- gd(X) = — 1/ log (rid^ 2 e 2 " 1 a 2 d ) , where A 
is an arbitrary dimensionless cutoff parameter, which is present due to the truncation 
in the perturbative expansion; it is the analogous of the cq parameter discussed in 
Subsection 13. II In the following is the mass of the bosons. Up to the second order 
in the running coupling, one can express the energy density in the following way 



£ = 



g d (X)(log (7r 5d (A)) -logATT 2 



(3) 



It is evident from the above expression that fixing the scattering length which 
appears in the definition of <7d(A) is not sufficient for determining the energy density 
if one does not declare its choice for A, that is the form of the coupling constant. 
Notice again that the coefficient of the second order term does depend on the choice 
of A. A convenient choice for simplifying the expression is to set A = e~ 27 /7r 2 , so that 
we can introduce the coupling <jy = — 1/ log (rt^a 2 ,) and we obtain 



2irh 2 n 2 d 
m d 



9d 



l + .9rf(log (vr.g d ) + 27+i 



(4) 



Now let us consider the case when the bosons are dimers, consisting of two paired 
fermions with mass m — rrid/2 and particle density n = 2n d - There must exist a 
regime where the binding is so tight that the equation of state of such composite 
bosons is the same as in the case of point-like bosons, with the simple replacement 
ad — > aa2D- Let us therefore introduce 77 = log {Uf^d) in the previous expres- 
sion, so that the composite bosons coupling turns to be gd = —l/\og(na 2 a\ D /2) = 
l/(log47r — 2rj — 2 log a). In such a situation the energy per fermion can be written 
in the following way: 



E 



2 



£p 1 

~2 2 gd 



l + .9d(log {irgd) + 27 + ^ 



(5) 
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(a) (b) 

Fig. 4. (a) g^i correlation function for the attractive gas at rj = 2.15 (g — — 0.25), 
compared to the short-range expression (b) g^i correlation function for the repulsive 
HD gas at g = 0.28. 



where — |£b|/2 is the contribution from the binding energy of the dimers. From a 
fit of this expression we obtain a ( i = 0.55(4)a2D, hi agreement with the four-body 
calculation in Ref. [12]. In order to exemplify the importance of the second order 
expansion ([5]), in Fig. [3] we also show the leading linear contribution alone, with the 
same choice of the coupling constant. Without the second order term it would be 
impossible to determine ad with accuracy. 

To our knowledge a theoretical analysis for the running coupling of a 2D fcrmionic 
non-relativistic fluid, analogous to the one in Beane Q]5] for the bosonic counterpart, 
is still lacking. It would be highly useful for interpreting the results of Subsection 
13.11 putting the comparison of the Quantum Monte Carlo data and the Fermi liquid 
equation of state on firmer grounds. 

The solid curve in Fig. [3] is a global fit of the data; the function to be fitted is 
a piece-wise defined function of 77, the matching point being a fitting parameter and 
the two matched functions being the ratio of second order polynomials of rj such that 
the global function is continuous at the matching point up to the second derivative 
and the extreme regimes coincide with the known perturbative results. 



4 Contact parameter 



The Contact parameter C is a property of short-range interacting gases, measuring 
the amount of pairing between the particles and relating a large number of observables 
with each other [20] . For example it can be obtained from the derivative of the equa- 
tion of state with respect to the coupling parameter C — (2irm/h 2 )d(nE/N)/d(log kp(i2D) 
or from the short-range behavior of the antiparallel pair distribution function g+i (r ) — > 

iC/kf, log 2 (r/a2_o) [H] (where r — > means R < r <^ l/kp, R being the range of 
the potential). 

We calculate the g-j-j. correlation function along the crossover using FN-DMC and 
extract the Contact by means of the following more refined formula, which better 
describes the behavior at small r (see Fig. 4(a) ): 



< ^ 4C 

9M r^ 4 



4i) 1 




(0) 



The additional terms permit a more precise determination of the Contact in the deep 
BEC regime, where the peak in at r — > is very narrow since a^D is quite small 
with respect to the mean interparticle distance. 
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Fig. 5. Contact parameter in the weakly repulsive regime. The solid line corresponds to the 
calculation from the derivative of the equation of state The triangles correspond to the 
extraction from the short-range behavior of g^i for the HD gas. 



The results for the BCS-BEC crossover have already been presented in Fig. 4 
of Ref. [5], where they have also been compared to the Contact extracted form the 
derivative of the global fit to the energy. The overall agreement between the two deter- 
minations of C is a useful check of the accuracy of the trial wavefunctions used in the 
FN-DMC approach. Small deviations in the region log(fcpa2D) ~ 1 point out the need 
of a better optimization of the trial wavefunctions, indicating that in the resonance 
region the JS wavefunction possesses too little pairing while the BCS wavefunction 
too much. 

In Fig. |4(b)| we also demonstrate the determination of C from the g^i of the 
repulsive HD gas. The short-range details are model dependent, nevertheless it is still 
possible to find an intermediate region which is universal. In this case the leading order 
formula for is used, since eq. © is valid only for r <C a,2D, which corresponds to 
the non universal region for the repulsive gas. In Fig. [5] we compare C extracted from 
g-fl to its analytical expression obtained from eq. ^fy: C/k F = [1 + (3 + 41og2)g]g 2 . 
Similarly to Fig. [5] deviations appear around g = 0.2. 



5 Conclusions 

To conclude we have detailed the functional forms needed for accurately extract useful 
information from the FN-DMC simulations of 2D Fermi gases. Both for the energy 
in the weakly interacting case and in the BEC regime and for the g^i correlation 
function the knowledge of the next-to-leading order correction is crucial in order to 
avoid ambiguity in the measured properties. We have also shown new results on the 
equation of state and the Contact of the repulsive gas in the weakly interacting regime. 

I want to acknowledge the support of the University of Trento and the INO-CNR BEC 
Center in Trento, where a large part of the results have been obtained, and the valuable 
discussions with S. Giorgini. 
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